Non-Markovian non-stationary completely positive open quantum system dynamics 
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By modeling the interaction of a system with an environment through a renewal approach, we 
demonstrate that completely positive non-Markovian dynamics may develop some unexplored non- 
standard statistical properties. The renewal approach is defined by a set of disruptive events, 
consisting in the action of a completely positive superoperator over the system density matrix. The 
random time intervals between events are described by an arbitrary waiting-time distribution. We 
show that, in contrast to the Markovian case, if one performs a system-preparation (measurement) 
at an arbitrary time, the subsequent evolution of the density matrix evolution is modified. The 
non-stationary character refers to the absence of an asymptotic master equation even when the 
preparation is performed at arbitrary long times. In spite of this property, we demonstrate that op- 
erator expectation values and operators correlations have the same dynamical structure, establishing 
the validity of a non-stationary quantum regression hypothesis. The non-stationary property of the 
dynamic is also analyzed through the response of the system to an external weak perturbation. 

PACS numbers: 03.65.Yz, 42.50.Lc, 03.65.Ta, 05.40.-a 



I. INTRODUCTION 

The theory of Markovian open quantum systems fH is 
well established from both a mathematical and a physical 
point of view. The theory of quantum dynamical semi- 
groups, casting the structure of completely positive (CP) 
trace-preserving maps, establishes that the Kossakowski- 
Lindblad equations are the most general admissible forms 
of evolution of the system density matrix. The applica- 
tion of these equations ranges from quantum optics 2] to 
quantum information theory . 

As far as the quantum non-Markov case is concerned, 
there exist different physical situations, and as a con- 
sequence a large variety of formalisms, from which a 
solid proposal for an approach to non-Markovian quan- 
tum dynamics [ll, 0] may emerge. A promising di- 
rection is afforded by the non-Markovian generaliza- 
tion of the Kossakowski-Lindblad equations. In a re- 
cent contribution Q , Barnett and Stenholm showed that 
the adoption of a time convolution between a mem- 
ory kernel and a Kossakowski-Lindblad operator, al- 
though appealing, may lead to unphysical results. How- 
ever, their attempt attracted the attention of many re- 
searchers to the search of the proper memory kernel 
for the time convoluted Kossakowski-Lindblad equations 
[i, 0, i, 1 M, EH, El El, El, El, El, E3 ■ The main focus 
of most of these papers has been devoted to the search 
of memory kernels that guarantees the CP condition of 
the solution map 0, 0, H i, El, El El, El, El • On the 
other hand, different microscopic interactions that lead 
to the convolution structure were established [H, 
and applied in the characterization of spin environments 
[13], quantum Boltzman equations with internal degrees 
of freedom [ll|, mesoscopic systems [1^, as well as to 



fluorescent systems coupled to complex self-fluctuating 
environments [2l| . 

It is worth pointing out that the adoption of a renewal 
approach based on the extension to quantum mechanics 
of the celebrated continuous time random walk (2^ leads 
naturally to the time convoluted structure that has been 
originally hypothesized by the authors of Ref. [1] , with no 
risk of violating the CP condition Q. In fact, the result is 
obtained through an average over infinitely many trajec- 
tories, each of which consists in a series of sporadic and 
consecutive transformations (system-environment colli- 
sions) of the system density matrix. The CP condition is 
trivially satisfied when each collision is written in terms 
of a CP transformation. 

The continuous time random walk formalism has be- 
come a fundamental tool when describing classical (non- 
Markovian) complex systems. In particular, the exis- 
tence of processes without a characteristic time scale 
(i.e., characterized by power law behaviors) has lead to 
an intensive review of the formalism and of its possible 
extensions. One of the recent motivations for studying 
that regime comes from the emergence of non-stationary 
phenomena in the fluorescence intensity produced by 
(blinking) nanocrrotal quantum dots under laser radi- 
ation [H, m, m, [2II . These experiments, have led many 
researchers to revisit some basic tenets and tools of equi- 
librium and non-equilibrium statistical mechanics, which 
are closely related the one to the other: the Onsager prin- 



master equations 
,[3l[33,[3a,[3E[33. 



29|, 



ciple [27|, [28!, |29| , non-stationary 
linear response theory [sl, ^1 [32 

Wiener-Kinchine theorem ^39] and the ergodic hypothe- 
sis 40.] • The extension of these principles and theoretical 
tools for dealing with non-stationary phenomena would 
be an outstanding breakthrough in statistical mechan- 
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ics. Besides the theoretical interest, there is an increas- 
ing number of experimental situations (4l| that can draw 
benefits from that theoretical progress. While there exist 
different issues that remain open, we are naturally chal- 
lenged to find a proper generalization to open quantum 
systems dynamics. 

The main goal of this paper is to show that non- 
Markovian CP master equations may fit non-standard 
non-stationary statistical phenomena, and then to ana- 
lyze the validity or extension of two cornerstones of the 
theory of open quanturn systems, i.e., the quantum re- 
gression hypothesis [3, H^j and linear response theory 
[set [Slj . The underlying dynamics of the system are de- 
fined by the renewal approach introduced in Ref. Q. 

To address the definition (in the context of an open 
quantum system theory) of a non-stationary decay (or a 
non-stationary quantum master equation), we introduce 
a system preparation at an arbitrary time posterior to 
the initial coupling between the system and the bath. 
With preparation we means an instantaneous CP opera- 
tion (like a measurement or any sudden CP transforma- 
tion) that leaves the system in an arbitrary state. Then, 
two times are introduced. One of them (t) measures 
the time at which the preparation occurs, and the second 
one (r) measures the time since the preparation. We shall 
use the term stationary decay (stationary master equa- 
tion) to denote a relaxation after the preparation done 
at time t that is independent of t, i.e., it only depends 
on T. On the same token, we shall use the term non- 
stationary decay (non-stationary master equation) to de- 
note relaxation processes whose form depends on <, i.e., 
its functional dependence on t is parametrized by t. In 
the Markovian case, the preparation always leads to the 
same stationary master equation. In the non-Markovian 
case, we show that even when the preparation is per- 
formed at arbitrary long times, the ensuing relaxation 
may or not reach a stationary regime. 

We also show that even in the presence of non- 
stationary effects, operators expectation values and cor- 
relations have the same dynamical structure, providing a 
generalization of the standard quantum regression theo- 
rem to a class of non-stationary quantum dynamics. 
The response of the system to an external weak pertur- 
bation, while it can be defined in terms of operator corre- 
lations [13, [Ml , generates strong deviations with respect 
to the Markovian case. 

The outline of this paper is as follows. In Sec. II, 
we review the renewal approach and show how the non- 
stationary effects arise. In Sec. Ill we obtain the evolu- 
tion of both operator expectation values and correlations, 
which allows us to establish a quantum regression theo- 
rem. In Sec. IV, the response to external perturbations 
is studied. In Sec. V we provide the conclusions. 



II. NON-STATIONARY DENSITY MATRIX 
EVOLUTION 

In the quantum application of the renewal approach 
the density matrix Ps{t) of an open quantum system 
S is determined by means of an average over an ensem- 
ble of infinitely many stochastic realizations, Pgir) = 
(/9^((t)), where (•••) denotes the average over realiza- 
tions and Pst(T) is the stochastic state associated to each 
trajectory. They consist of a sequence of disruptive (col- 
lisional) events occurring at random times. The times 
elapsed between two consecutive events are randomly 
drawn from a waiting-time distribution density w{t), sat- 
isfying w{t) > 0, and /p w{t)dt — 1. Each event is asso- 
ciated to an arbitrary CP transformation £ of the system 
state. It is defined by the Krauss form 

£[p]^Y.^.^P^l (1) 

i 

where p is the system state prior to a given event. The 
operators Ci satisfies the condition C|Ci =1. Further- 
more, here we assume that between consecutive events, 
the evolution of the system is defined by the propagator 
exp[f>C5]. The superoperator £5 is the Liouville super- 
operator and corresponds to a unitary transformation. 
Nevertheless, we remark that most of the results hereby 
developed also apply when £s is a standard Lindblad 
superoperator, i.e., when the evolution between events 
corresponds to a Markovian (CP) dissipative evolution. 

By construction, each realization, and as a consequence 
the average over the realizations, guarantees the CP con- 
dition of the solution map. The system's dynamics begin 
at time t = Q (system-environment coupling). As stressed 
in the introduction, we let the system evolve up to the 
time t > that we set to be the new origin of time. 
The earlier work of Ref. [3| is confined to the condition 
t = and the main aim of this paper is to solve the non- 
stationary issues raised by the condition i > 0. The main 
idea of the method that we use is as follows. First of all 
we study the time evolution of pg from ^^(0) to ^^(t-l-r), 
and we interpret the exact expression of pg[t -\- r) as the 
the density matrix Ps{t) that will be expressed in terms 
of the initial condition p{t). It is straightforward to get 
the exact expression of Ps{t)^ which reads 

00 pT+t 

Ps{r) = J2 dt'Vo{r + t~t')W^-Ht')ps{0). (2) 

The superoperator >V(")(t) is defined in the Laplace do- 
main, (t u, t —^ z) as 

W^''\u) = [£w{u-/:s)r, (3) 
while the superoperator Vq{t) reads 

Vo{u) = Poiu^Cs), (4) 
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where Po{u) is the survival probability associated to 
'w{u), i.e., 



Poiu) ^ 



(5) 



The expression given by Eq. ([21) is a sum over all possible 
realizations, each of them corresponding to a stochastic 
process with n collisions. We have now to express it in 
terms of the initial condition 



By using recursively the relation 

n;(")(r + t) ^ r dr'Wir - t')W^''-^\t' + t) (7) 



Jo 



we rewrite Eq. ([2]) as 

Ps{t) = So(r,t) + E / dr' / dr" (8) 
xPo(t - t')>V("-i)(t' - r")S^(T", t), 



where we have defined 



So(T,t) = V / dt'Vo{T + t-t')W^-\t')ps{0), (9) 

oo „( 

S^(T,i) = V / di'W(T + t-t')W(")(i')Ps(0). (10) 
The Laplace transform of pg{T) of Eq. ([8]) then reads 

oo 

Ps(u) =So(7/,z) + E7'o(u)W^"~'^(w)S„(u,z). (11) 

n=l 

By using the relation 



/ dr / (ite-"^e~^*(7(T + 1) 
Jo Jo 



g(u) - g{z) 



, (12) 



which is valid for any arbitrary function g{t), with the no- 
tation g{u) = Jp°° dTe~"'^ g{T) and g{z) = dte~^*g{t), 
the Laplace transforms of So(r, t) and Si„(t, i) are writ- 
ten as 

„ , , >V(u)->V(z) 1 , , 

Here, Ps{z) is the Laplace transform of Ps{t), given by 
Eq. ([6]). By plugging these expressions into Eq. pT|) . 
after some algebra we get 



Ps{u) = Q{u)[ps{z) + Ip(u, z)]. 



(15) 



where the propagator Q iu) is defined by 



(16) 



u-[Ls+ LK(u - Cs)^ 
and the inhomogeneous term by 

Ip(u, z) = C^{u ~Cs,z- Cs){z ~ Ls)Ps{z). (17) 
In the time domain, Eq. (|15p becomes 

rfp5(T) 



i^sPsir) + / Vi^(r - T')/:e(^--')^-ps(T') 

"'0 



+lp{T,t). 



(18) 



This equation is one of the central results of this section. 
It defines the evolution of the average density matrix of 
the system in the interval (t, t + t), with the initial con- 
dition Eq. ((6|). The superoperator C is defined by 

C = £-l, (19) 

which in turn can be written with the Lindblad structure 

m^\Y.{[C,,.C}] + [C,.,Cl]}. (20) 

i 

The memory kernel function K{t) is defined in the 
Laplace domain by 



K{u) = 



uwiu) 
1 — w{u) 



(21) 



The inhomogeneous contribution lp{u,z) [Eq. (fTT]) ] is 
proportional to the function 



A(u, z) 



w[u^z) 'w{u)/z 



1 — w{u) 1 — w{u) ' 
where the function z) reads 

w{u) — w{z) 1 



w{u, z) 



z — u 1 — 'w{z) 



(22) 



(23) 



By using the relation of Eq. (|12p . w(m, z) can be written 
in the time domain as 

oo „f 

u;(t,<) = w(r + i) + V / dt'w{T + t-t')w^''\t'), (24) 

where w^'^\z) = [w(z)]". This expression allows us to in- 
terpret w{T,t) as a conditional waiting-time distribution 
density, and more precisely as the probability distribu- 
tion density of meeting the first event at time t, given 
that the observation time (of events) begins at time t. 
The second term on the right hand side of Eq. takes 
into account all possible events at times earlier than t. 
Consistently, notice that for t = 0, w{t, 0) = w{t). 
In the time domain, the function A(m, z) read 



A(r,t) = /(r,i)-/(r,0), 



(25) 
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where we have introduced the (sprinkhng) distribution 

w{u, z) 



1 — w(u) 

By writing this expression in the time domain, 



(26) 



n=l 



dT'w'^'^\T^T')w{T\t), (27) 



it follows that /(r, t)dT is the probability of an event oc- 
currence in the time interval (r, r + dr), given that the 
observation time begins at time i, regardless of whether 
or not any event occurred at earlier times. It satis- 
fies the relation /(t, 0) = w(0,t). On the other hand, 
notice that the function T(m, z) = zA(u, z) [appear- 
ing in Eq. (|17p ] can be written in the time domain as 

T{T,t) = {d/dt)f{T,t). 

Both the kernel K{t) and the inhomogeneous contri- 
bution Ip(r,t) are clear signatures of the non-Markovian 
property of the evolution Eq. ^5]) . Using Eq. p5|) . it 
is easy to realize that the evolution of Ps{t) can always 
be rewritten as an homogeneous evolution [see for exam- 
ple Eq. ([3T|) ]. Nevertheless, the inhomogeneous structure 
allows us to understand which is the effect of shifting 
the initial time condition from ^5(0) to Ps{i)- In fact, 
Eq. (fT7)) tells us that in the interval (i, t+r) the departure 
of the system time evolution from its time evolution in 
(0, t) is measured by Eq. ([25]) . Consistently, for t = 0, the 
inhomogeneous contribution vanishes, i.e., Ip(T,0) — 0. 

No departure of the system evolution in (t, t + t) 
from the time evolution in (0, t) must occur in the 
Markovian case. This case is recovered by assuming 
an exponential waiting-time distribution (Poisson case), 
w{t) = 7exp[— 7r]. From Eqs. (|23p and (|26|) it follows 
that w{T,t) — w{t) and f{T,t) — 7, thereby implying 
the vanishing of the inhomogeneous term, and the rela- 
tion K{t) — 7(5(t), which turns Eq. pB]) into a standard 
Lindblad equation 



dr 



[Cs + jC]psiT). 



(28) 



On the other hand, we remark that in Eq. psp . and in 
Eq. (|17|) as well, the superoperator C may be assigned 
an arbitrary Lindblad form. In fact, when C 7^ £—1 [see 
Eq. p^ ] the superoperator £ [Eq. ^] can be defined 
as £[p] = {I -I- [e"'^ — I]}/9, where k must be interpreted 
as a control parameter. Then, we recover Eq. ITSl) with 
an arbitrary C in the limiting condition in which simul- 
taneously K — > and the number of events per unit of 
time go to infinity, the last limit being controlled by the 
distribution /(t,0) of Eq. 



A. Initial preparation at time t 

The initial condition associated to Eq. ^5)) is given by 
Ps(t), Eq- ©J which in turn carries information about 



the system dynamics in the interval {0,t). Therefore, 
Eq. (jlSp does not give more information than a master 
equation describing the evolution in the interval (0, t+r). 
Nevertheless, the master equation ([T5|) may acquire a dif- 
ferent status if the initial condition at time t can be cho- 
sen as any non-equilibrium form of the density matrix 
Pg. This is done by introducing the main ingredient of 
our formalism, i.e., by adopting the concept of prepara- 
tion, namely, a change pg{t) ^ pjj, compatible with a 
CP transformation 11. The role of the preparation is to 
erase the dependence of the evolution on the previous 
history of the system without erasing the memory of the 
universe, i.e., the system-environment arrangement. 
In the Laplace domain the preparation is defined by 



(29) 



which in the time domain yields pg{t) — > pYiexp[Cst\. 
The extra unitary contribution is introduced to take into 
account that the Hamiltonian evolution defined by £s 
begins at time t = 0. Similarly, one can interpret Eq. ([^5]) 
as a preparation in an interaction representation with 
respect to As- 
under the preparation condition of Eq. (|29|) , the time 
evolution structure of Eq. psp is still valid, provided that 
the initial condition is fixed to be pjjexp[£sf], with the 
inhomogeneous term now reading 



Ip(r, t) = £A{t, t) exp[(T + t)£s]Pn- 



(30) 



We remark that the time evolution of Eq. ifTSl) with the 
contribution of Eq. (f30|) remains a CP structure. In fact, 
also its solution admits an interpretation in terms of tra- 
jectories that preserve the CP condition. 

The time evolution of the density matrix generated 
by Eq. with the inhomogeneous contribution of 

Eq. (j30[) . in principle depends on the preparation time 
t. Which is the form of the dependence of the density 
time evolution in (t, t + r) on t? Of particular interest is 
to assess under which conditions this dependence is lost, 
so as to generate in the long-time limit the stationary 
behavior defined in Sec. [H In the case where the solu- 
tion becomes asymptotically stationary, it is of interest to 
assess if this stationary time evolution (i — > cxd) is charac- 
terized by non-Markovian effects stronger or weaker than 
the time evolution with the preparation stage coinciding 
with the initialization stage, i.e., with t = 0. These im- 
portant questions will be answered with the help of the 
simple examples discussed in Section [II Bl 

As a last but not least remark of this Section, let us 
notice that the time evolution of Eq. , with the prepa- 
ration condition of Eq. (j29p , can be easily written in an 
equivalent form, as an homogeneous time evolution, as 
follows 



dps{T) 
dr 



^sPsi'r) 



dr' 



dT"Mt{T-T') 



Kt{T'-T")Ce^^'-^"^^^Ps{T"), 



(31) 
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where A^2('u) = [1 + CA{u- Cs, z- Cs)]~'^ , and Ktiu) = 
uw{u, t)/[l — w{u)]. While this expression avoids the com- 
phcation arising from the presence of an inhomogeneous 
term, the kernel structure is more complicated, involv- 
ing all powers of the Lindblad superoperator C. Eq. ((3T|) 
recovers and generalize the classical master equation ob- 
tained in Ref. (H. 



B. Examples 

To make more transparent the spirit of the renewal ap- 
proach of this paper, here we illustrate it in action on two 
exemplary cases of the same simple model. While the mi- 
croscopic origin the superoperator £ and the waiting time 
distribution w{t) is not completely understood 0, IT5| . 
from Eq. (fT5|) it becomes clear that the former object 
defines the underlying Lindblad like structure Eq. (pn|) . 
Then, it establishes the coupling between the density ma- 
trix elements. On the other hand, wit) can be settled in a 
phenomenological way as a function of the characteristic 
system decay behavior (see next examples). It leading 
property is the average waiting time, t'w{t)dt, which 
may be finite or divergent, the last case giving rise to 
strong non-stationary effects. 

As a simple model, we consider a degenerate two-level 
system (£5 0) and the superoperator 



(32) 



The time evolution of the expectation values of the Pauli 
matrixes cr^, Siir) = TTs[ps{T)ai], with i = x,y,z, is 
given by Eqs. and (001) [or equivalcntly by Eq. CT ]. 
and it reads 



dr 



dT'kt{T-T')Sx.Y{T'), (33) 



while Sz{t) = 5*^(0). With S'i(O) we denote the expecta- 
tion values after the preparation. The kernel is defined 
by its Laplace transform Kt{u) = uw{u,t) /[I — w{u^t)\. 
The solution of Eq. ^ is 



(34) 



Sx.Y{T) = Sx.Y{Q)Po{T,t). 



where PQ{T,t) is the survival probability associated to 
w{t^ t), i.e., Po{t, t) = 1 — /g^ dT'w{T' , t). It can be rewrit- 
ten as 

Po{T,t)^Po{T + t)+ f dt'Po{T + t-t')f{t',0). (35) 



Here, Pq{t) is defined by its Laplace transform Eq. (O, 
while f{t,0) follows from Eq. (EH), i.e., /(z,0) = 
w{z)/[l — w{z)]. 

The decay of the expectation values of Eq. ([M)) de- 
pends on both T and t. Its explicit analytical form de- 
pends on the choice done for the waiting-time distribu- 
tion 'w{t). As first case, we select the bi-exponential case 



w 



(t) = Palae-''-' + Pblb^-'''\ 



with Pa + Pb = l- The distribution Eq. ^ {t = Q) reads 
/(r, 0) = (7) e{T) - [(7> - {r)-\l e-^n, (37) 

where 9{t) is the step function and we have introduced 
the parameters (7) = Pa7a + Pblb^ (t) = Pala^ + 
Pblb^ = /o°° Tw{T)dT < 00, and 77 = Pa7h + Pbla- No- 
tice that after a transient of order I/77, the sprinkling 
distribution, as in the Markovian case, is constant, i.e., 
/(t,0)~1/(t) >0. 

The coherence decay, independently of the time t, can 
be written as 



PoiT,t) ^ Pa{t)e-''^^ + Pb{t)e-''^\ 



(38) 



All the dependence on the preparation time is carried out 
by the weights Pa{t) and Pb{t). Their explicit form fol- 
lows straightforwardly from Eq. ([55)1 as a superposition 
of exponential functions. They satisfy the boundary con- 
ditions Pa{0) = Pa, and Pb{0) = Pb- In the limit i — * 00, 
the asymptotic stationary decay reads 



Po{t,oo) = 



Pa 



Pb 



(t) lb 



(39) 



In Fig. 1, we show the decay defined by Po{T,t) for dif- 
ferent preparation times t, and for two different sets of 
characteristic parameter values. After a transient of or- 
der ?7, both cases reach a stationary decay regime. By 
comparing these figures one with the other, we realize 
that the asymptotic decay may yield arbitrary depar- 
tures from the dynamics generated by setting t — 0. In 
fact, in Fig. la, the asymptotic decay is almost expo- 
nential while the initial one is bi-exponential. In Fig. lb 
the inverse situation is observed. This simple example 
demonstrates that no general conclusion can be drawn 
about the properties of the stationary time evolution. 

As a second exemplary case, we consider the waiting- 
time distribution 



w{u) 



(40) 



where the units of Aa are 1/sec", and < a < 1. Note 
that for a = I this expression reduces to the Laplace 
transform of an exponential function. The kernel Eq. (|21D 
read K [u) = AaU^~°'. As is well known 44], this kind 
of kernel is related to a fractional derivative operator. In 
contrast to Eq. ([57)) . here we get 



/(^,0) 



An 



r(a) t(i-«)' 



(41) 



(36) 



where r(a;) is the Gamma function. Then, in this case, 
lim^_,oo /(■'■, 0) = 0. This property is directly related to 
the divergence of the average period between events, i.e., 
Tw{T)dT = cx) [4i|. 
By using Eq. ([55)1 and the fact that ^i(''') is, in 
this case, a Mittag-Leffler function 0, Ej], we can 
write PQ(T,t) [45] under the form of a series expansion. 
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FIG. 1: Coherence decay Po{t, t) [Eq. p4|) ]. for different times 
t, for the waiting time distribution Eq. (|36p . In (a) the param- 
eters are Pa = 0.8, Pt = 0.2. In (b) are Pa = 0.99, Pt = 10"^ 
In both cases the rates are 7^ = 1, 75 = 0.05. Both, r and t 
are measured in arbitrary units (a.u.). 



By using the property that for AaT°' ^ 1, Po{t) ~ 
Aa/[T°'T{l — a)], when t ':$> twe get the following asymp- 
totic expression 



r(i - a) 



(r -I- t)" aT{a) (r + i)" 



(42) 



Therefore, in this case there not exists an asymptotic 
stationary decay. In fact, this expression shows that at 
any time the decay dynamics depends on the preparation 
time t. In Fig. 2 we plot the function Pvi{t, t), Eq. (|35p . for 
different values of t. Eq. (|^^ correctly fits their asymp- 
totic decay behavior. Consistently, we found that there 
not exists a stationary decay behavior. 



III. REGRESSION HYPOTHESIS 



The generalization of the classical regression hypothe- 



sis 



2l\ to a quantum context is called quantum regression 
theorem 0, \^ . It states that operator expectation val- 
ues and operator correlations have the same dynamical 
behavior. Here, we explore the possibility of generalizing 
this theorem to the renewal case. 



10° d 
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FIG. 2: Coherence decay Poir, t) [Eq. pi)) ]. for different times 
t, for the fractional waiting time distribution Eq. (|40p . The 
parameters are = 1, a = 1/2. Both, t and t are measured 
in arbitrary units (a.u.). 



A. Operators dual evolution 

In order to define operator correlations, we have to 
move from the Schrodinger to the dual or Heisenberg rep- 
resentation. In the renewal case here under study this 
corresponds to convert the stochastic time evolution of 
the density matrix Pgt(t) into the stochastic time evo- 
lution of operators. All this rests on the fundamental 
relation 



Ait) = Tvs[A{0)ps{t)] = Trs[ps{0)A{t)], 



(43) 



where the mean value A{t) of a system operator A can be 
written in terms of the initial density matrix Ps{0) and 
of the evolved operator A{t). 

Let us define the dual superoperators and by 
the relations 



Tr5[Ae*^^p] 
TrsiASp] 



Tr s[pe*^s A]. 
Tr s[p£* A]. 



(44a) 
(44b) 



Eq. dH) yields the (averaged over realizations) operator 
time evolution 



A{t 



00 pT + t 

,1=0 -^0 



dt'w*'^''\t')r^{T + 1 - t')A{o). 

(45) 



Here, W*'^"\z) = [w{z - C*)£*Y' arises from Eq. ©, 
and 'Pf{z) = Po{z - C^) from Eq. Note that, as 
a consequence of the prescription of Eq. ((43)) the W'^ 



superoperator applies after the one rather than before 
it as in Eq. 

In conclusion the dynamics of the stochastic operator 
Ast resembles that of the stochastic density matrix p^^. 
It consists, too, of time intervals with the time evolu- 
tion driven by exp[t£^], and of others, corresponding to 
the action of the superoperator E'^ , where it is driven 
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by disruptive events. Nevertheless, notice that when 
[Cg.E'^] ^ 0, the time ordering of the superoperators 
is reversed as a consequence of turning the Schrodinger's 
into the Heisenberg's representation. 



B. Operator expectation values and correlations 

For the main purpose of working with simplified ex- 
pressions, in this section we make all calculations in the 
interaction representation with respect to £g, and we 
assume that 

[4,f#]=0, (46) 

which in turn in the Schrodinger representation yields 
[Cs,£] = 0, or equivalently the commutation condition 
[CstC] = 0. As a consequence of this condition, Eq. ((45|) 
(Eq. (I2|) makes the operator (density matrix) stochastic 
dynamics, in the interaction representation with respect 
to (Cs), only consist of the application of the super- 
operator {£). Notice that all the expressions obtained 
in the previous section, in the interaction representation, 
remain valid by setting Cs 0. 

The operator expectation values, in the interval (i, t + 



t) are written as [A(t + t)^ A(t)] 

A(7) = Trs[ps{0)A{T + t)], (47a) 
= TTs[ps{0)CiA{r,t). (47b) 

With A — {Ai, A2,- ■ ■)'^ , we denote a vector of system 
operators defining a complete basis in the dual (opera- 
tors) space. I denotes the system identity operator. The 
operator correlations are written as 



0(i)A(t + T) = Tvs[ps{0)O{t)A{t + T)], (48a) 
= TTs[psiO)CoAiT,t)]. (48b) 

Here, O denotes an arbitrary system operator. The aux- 
iliary operators function CuviT,t), acting on arbitrary 
system operators U and V, are defined by the expression 



Cuv{T,t) 



OC oc 

EE 

m— n— 



P{T,m;t,n){£*)"[U{£*)'"'[V] 



P(r, to; i, n) is the probability that n events occur in the 
interval (0, t) and to events in the interval {t, t+r). Under 
the condition (|46p . in the interaction representation with 
respect to Eq. (|49| follows straightforwardly from 
the stochastic dynamics associated to Eq. (|45| . 

The set of probabilities P(t, to; t, n) can be written as 



dt'Pa{T + t-t')w^''\t'), (50) 



P(r,0;t,n) = 
when m = 0, and as 



P(T,TO;<,n) = / dT'Pa{T-T') dr" w^"'-^\t' - t") 



for m > 1. As before, the function M;(")(r) is defined by 
its Laplace transform w'^"^(u) = [w{u)]'^ . In the Laplace 
domain, {t u, t ^ z), after using Eq. (fT2|) . we get 

„/ N 1 — zw(u, z) „ , , r / s-,^ 

u 

and for to > 1, 

P{u,m; z,n) = Pq{u)[w{u)]'"~^ zw{u, z)Po{z)[w(z)]" , 

where w{u, z) is defined by Eq. ([25)) . 

After some algebra based on Eq. we write 

Cuviu, z) = g*{z)Ug*{u)[l + zA{u, z)£*]V. (52) 

The function A(it, z) is defined by Eq. ^ and g*{u) 
denotes the propagator 

1 



g*{u) 



(53) 



-K{u)C#' 

where = £'^ — 1 is the dual superoperator asso- 
ciated to £, Eq. (fT9)) . The kernel K{u) follows from 
Eq. (HJ). Therefore, taking into account that these ex- 
pressions were derived in an interaction representation 
with respect to Cs, we obtain that g'^{u) is the dual 
propagator associated to g{u), Eq. p^ . 

From Eq. (|52p . after introducing the density matrix 

Ps{z) = g{z)ps{G), (54) 

the mean values Eq. (|¥71) and correlations Eq. (HS)) read 

A(7) ^ TYs{ps{z)g*{u)[l + z/^{u,z)C*]K}, 

0{t)A{t + T) ^ i:Ys{Ps{z)Og*{u)[l + zA{u,z)C*]A}. 

Here, for the sake of shortening the notation, we use the 
symbol = to indicate that the left and right hand side of 
the equality are written in the time and Laplace domain, 
respectively. These equations yield the desired expres- 
sions for operator expectation values and correlations. 
They can be straightforwardly written in terms of den- 
sity matrix propagators as 



A(r) 



Tis{Ag{u)[l 

i:rs{Ag{u)[i 



0{t)A{t + T) 

In the Markov case, i.e., when K(u) 



C/\{u,z)z]ps{z)}, (55) 
C/^{u, z)z]ps{z)0}.{bQ) 
7, these results 



recover the expressions that follows from a microscopic 
derivation based on a Born-Markov approximation 
Furthermore, by using the same calculations steps it is 
possible to demonstrate that 



X / dt'w{T" + t-t')w^^''\t'), (51) 
Jo 



0{t)A{t + T)0{t)^Trs{Ag{u)[l+zA{u, z)C]Opsiz)0}. 

(57) 

From Eqs. (f55|) and ([56|) it is immediate to realize that 
expectation values and correlations have the same dy- 
namical structure, showing that the classical Onsager re- 
gression hypothesis can be extended to this context. We 
make this fact even clearer by introducing a preparation 
at time t [Eq. thereby implying the transformation 

zpg{z) — > Then, the preparation can be interpreted 
as a sudden fluctuation at time t. The earlier expressions 
indicate that the operator correlation dynamics depend 
on the dynamical decay of this fluctuation. 
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C. Evolutions 

We can explicitly show that operator expectation val- 
ues and correlations have the same dynamical behavior. 
Here, we obtain the inhomogeneous equations of motion. 
Nevertheless, as in Eq. ipij) . they can be rewritten as ho- 
mogeneous ones. By defining a matrix M by the relation 



Ti-s[ACO] = MTr5[AO] 



(58) 



which acts on the indexes of vector A, from Eq. ([55]) it 
is possible to get the evolution 

— A(r) = / dt'K{T - t')MA(t') + riA(T, t), (59) 
dr Jo 



while from Eq. ([56)l . for the correlations it follows 

— 0{t)A{t + T) = [ dT'K{T - T')MO{t)A{t + t') 
dr Jo 

+roA{T,t). (60) 

The inhomogeneous terms riA(T,t) and ToA{T,t), tak- 
ing into account the preparation [Eq. ((29l) ]. follow from 



ToAiT,t) = A(T,t)MTrs[OApn]. 



(61) 



These expressions show that even in the presence of 
strong non-Mar kovian non-stationary effects, the regres- 
sion hypothesis is still valid ^2^, 2M ■ 



Note that G{u) and lp{u, z) are defined by Eq. (|16p and 
(|17p . respectively. After performing the preparation at 
time i, the inhomogeneous term follows from Eq. (|30p . 

The work of Ref. discussed the validity of the re- 
gression hypothesis in the context of non-Markovian dy- 
namics based on Lindblad rate equations [lB| . It was also 
found (see Sec. 6) that an external (non-commuting) field 
breaks its applicability. Nevertheless, in contrast with 
the present formalism, the non-Markovian effects admit 
an underlying Markovian description. Furthermore, the 
regression theorem was studied by analyzing the dynam- 
ics at the initial and at the asymptotic time. In that case, 
the vanishing of the inhomogeneous term is a necessary 
condition for the validity of the regression hypothesis. In 
spite of these differences, both formalisms lead to consis- 
tent and non-contradictory results. 

We remark that similar conditions but no equivalent 
were found in different contexts. In Refs. 
|48| the validity of the quantum regression hy- 
pothesis beyond a weak coupling regime was discussed 
[T6| . In Ref. [i^, a commutation property between the 
system-bath interaction and the system operators was 
derived by using a stochastic wave vector formalism and 
taking into account a Bosonic bath described in a rotat- 
ing wave approximation. All these results suggest that, 
beyond a Markovian regime, the validity of the quantum 
regression hypothesis strongly may depends on the un- 
derlying microscopic dynamic. Nevertheless, the search- 
ing of general applicable criteria should not be discarded 

0. 



D. Discussion 



IV. LINEAR RESPONSE THEORY 



The previous analysis demonstrates that the condition 
(|46p guarantees the fulfillment of the quantum regression 
hypothesis. This constraint is satisfied, for example, by a 
two-level system with Liouvillian >Cs[»] — —ii^A[(^z, •]/2, 
where ctz is the z-Pauli matrix, lua its transition fre- 
quency, and £ (or equivalently£[»]) define a dispersive 
or thermal reservoir (see respectively Eqs. (71) and (74) 
in Ref. Q)- When the two-level system is subjected to 
an external field the regression hypothesis may be bro- 
ken. In fact, when the condition ((46|) is not fulfilled, 
the regression hypothesis does not hold true in any case. 
Nevertheless, by writing 



Cs = Co + eCi, 



(62) 



where [£o,^-] = 0, and 7^ 0, it is possible to prove 

that to first order in the parameter e, the regression hy- 
pothesis, independently of the specific structure of £1, is 
still valid. In the Schrodinger representation, the opera- 
tor expectation of Eq. (|55p and the operator correlation 
of (|56p. to first order in e, read 



A(r) « TTs{Agiu)[ps{z) + lp{u,z)]}, (63) 



Here we analyze the response of the system, whose den- 
sity matrix evolution is given by Eq. (jlSp . to an exter- 
nal time dependent perturbation. In the stationary case 
[30I [sH , the system response to weak external perturba- 
tions is expressed in terms of response functions that are 
proportional to the cross correlation function between the 
variable of interest and a system variable coupled to the 
external field. Here, we show that a similar result can 
be established, but that, nevertheless, strong departures 
from the predictions of the stationary theory may arise. 

To simplify the analysis, in the following calculations 
we assume that i = 0, i.e., that the system-environment 
coupling (initialization of the renewal dynamics) coin- 
cides with the preparation time, and that the coupling 
with the external field is switched on at the same time. 
In this case the absolute time coincides with the distance 
T from the system-environment coupling. The average 
system state Ps(t) is written as a series in the external 
perturbation 



(65) 



0{t)A{t + T) « TTs{Agiu)[ps{z) + lp{u,z)]0}. (64) 



The parameter A measures the strength of the external 
perturbation, p^g"^ (r) corresponds to the dynamics in the 
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absence of the perturbation. The contribution p's\t) 
can be obtained from an average of the perturbed real- 
izations or from the (perturbed) master equation defining 
the density matrix evolution, Eq. (jlSp . In general, one 
may assume that the external perturbation affects either 
the unitary or the dissipative dynamics. Hereby we ana- 
lyze both cases. 



A. Perturbing the dissipative dynamics 

Here, we consider the case when the external pertur- 
bation affects (or is coupled to) the dissipative dynamics. 
The two contributions in Eq. (j65p are evaluated by aver- 
aging the perturbed stochastic realizations. The zcro-th- 
order contribution reads 



Ps\r) = J2 / dr'7'o(r-r')W(")(r')P5(0). 



(66) 



Notice that this expression follows straightforwardly from 
Eq. ([2]) . The first-order contribution is determined by an 
average over all possible trajectories, in each of which the 
external perturbation acts only once, to fit the request of 
a linear response. Then, p^g\T) becomes the double sum 



Ps 



f'T />ri f'T2 

xWW(ti - T2)O(T2,r3)W(™)(T3)Ps(0)- (67) 



Each sum takes into account all possible events, pre- 
ceding and ensuing the action of the external perturba- 
tion. In this expression, Vq (t) and W(")(t) are defined 
by Eqs. ^ and ^ respectively. The influence of the 
external perturbation is described by the superoperator 
0{t2,T3). This is done either by making the superoper- 
ator £ change with time, without affecting the times of 
event occurrence, or by allowing the external stimulus to 
slightly change the times of event occurrence, namely a 
little bit earlier or later, according to the system's state. 



1. Perturbing the event superoperator 

In Ref. [sgJ the linear response theory was analyzed on 
the basis of a classical two-level system where the pertur- 
bation does not affect the time of occurrence of an event, 
but that the coin tossing selecting the fluctuations sign 
has a time-dependent bias. In Section |TT] we have seen 
that the Krauss operators £[p] signal the occurrence of 
renewal events. To realize a perturbation on the system 
of the same nature, we have to assume that the time of 
occurrence of collisional events is not affected by the ex- 
ternal perturbation, but that the specific form of £[p] is. 
Then, the Krauss operator is written as 



where the superoperator C(t) satisfies Tis{0{t)[p]} = 0. 
For simplicity, we assume 



0(r)=e(r)0. 



(69) 



^(r) is a scalar function that defines the temporal depen- 
dence of the external perturbation. The superoperator 
0{t2,T3) appearing in Eq. (|67[) can then be written as 

0(t2, Tg) = 0{t2) exp[(T2 - T3)£sMt2 ~ T3). (70) 

By working in the Laplace domain on the contributions 
to each sum of Eq. (I67|) . after some algebra, we get 



ps{t) ^ g{T)ps{o) + x dT'g{T-T')o{T') (71) 

Jo 

X r dT"K{T' - T")e'^^'^^"^^^g{T")psiO). 

The propagator Q(t) is defined in the Laplace domain by 
Eq. nil), while the kernel K{t) is defined by Eq. (|2T|l . 

Eq. (|7T|l generates the system's response to first order 
in the perturbation strength. A, and consequently the 
system's linear response. For the operator expectation 
values, after introducing the assumption ((69)) . we get 



A(t) - Ao(t) + A / dr'x^oir, r')e(T'), (72) 

JQ 



where the zero-th-order contribution reads Ao(t) — 
Tts[-^Q{t)Ps{^)]: and the response function is given by 

XAo(r, r') = Tr5[Ag(r - r')Op/(r')]. (73) 

With Pfir'), we denote 

Pf{r')^ r dr"K{r'~T")e(-'~-"^^^g{r")ps{0). (74) 



Evidently, the response function Xao(''"i ''"') ^as the 
structure of an operator correlation [see Eq. ((57)) ] . This 
fact becomes more evident when the initial density ma- 
trix corresponds to the stationary state p'^ of the unper- 
turbed evolution, i.e., 

^ hm (?(r)ps(0), (75) 

r— *oo 

and the condition 

Cs[p^] = Ps, (76) 
is satisfied. Then, the expectation values are written as 



£{T)[p]^£[p] + XO{T)[p], 



(68) 



A(r) =A^ + \ dr'xIo(r,T')e(r'), (77) 

where Aqo = Trs[Apg']. The response function becomes 
x'S:a{r, t') = Trs[Ag{T - r')Op^]/(r', 0). (78) 
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Here, the function /(r,0) is defined by Eq. ([271). The 
quantum statistical average over can be read as a 
correlation between the operator A and the superopera- 
tor O [see Eq. dST]) with zpg{z) p'f and A(m, z) 0]. 

For Markovian dynamics [Eq. (pS)) ]. where f{T',Q) 
7, the response Y'?'/n(T,T^) depends only on (r — r'). 
This stationary condition [30, |3l| is broken in the non- 
Markovian case. In fact, the presence of the factor 
/(t',0) implies that XAoi'^^'''') depends separately on 
both T and t'. Depending of the behavior of /(r',0) 
[see Eqs. (|37|) and (|4T|l ] in the long-time regime the sys- 
tem may become insensitive to the external perturbation. 
This effect, some times called death of linear response 
[ssj l , was found in classical systems in Refs. [sl, [s^l . The 
present analysis leads us to conclude that the same amaz- 
ing phenomenon may be observed in quantum systems. 

The previous results also follows from calculations 
based on the density matrix evolution, Eq. ((T8)l . On 
the other hand, we assumed that the external pertur- 
bation is switched on at the initial time. The general 
case, i.e., when the renewal dynamics start at time t = 
and the external perturbation is switched on at time 
t > 0, can be discussed following the same calculations 
steps. The final expressions involve some extra contri- 
butions. Nevertheless, under the assumption ((62)l . to 
first order in e, Eq. (|7ip remains valid under the re- 
placements Ps{Q) ^iid ^i''') ~^ ^t{T), where 
Kt{u) = uw{u,t)/[l — w{u)]. Similarly, Eq. (|75|) remains 
vahd under the replacement /(r', 0) /(t', t). 

As an example we consider a two-level system, with 
states {|±)}, whose unitary evolution is defined by 
'^S'M — ~i^[<^x, p\/'2, where ax is the x-Pauli matrix 
in the basis {|±)}. The perturbed superoperator read 

a, 6 

where (Tab = |a) {b\ , (a, b) = ±, and A|f (t)| < 1. It can be 
rewritten as £(r)[p] = [l+X^{T)az]/2. Then, both the un- 
perturbed and the perturbed terms turn out to be inde- 
pendent of the initial state p. The unperturbed dynamics 
correspond to a depolarizing channel 3] , with stationary 
state pg' — 1/2. The perturbed dynamics modulate the 
probability of transitions between the two states [IHl • 

From Eq. ([TZ)) . the difference between the upper and 
lower populations, i.e., the mean value of the z-Pauli ma- 
trix cr^, Sz{t) = TTs[pg{T)az], reads 

Sz{r) - A /V'Po(r - r') cos[r!(T - t')]/(t', O)e(T'). 
Jo 

(80) 

As the transformation £(r)[p] does not depend on p, it is 
simple to prove that Eq. ([80| also corresponds to the 
exact solution to all orders in A. On the other hand, 
notice that the results obtained in Ref. [1^, are re- 
covered in the limit 0. Using Eq. ([79)) . it is easy 
to establish the nature of the stochastic dynamics asso- 
ciated to Eq. dSni), i.e., Sz{t) = {Sstir)), where (• • • ) 




FIG. 3: Mean value of the z-PauIi matrix (full line), Eq. JHOj, 
driven by a periodic external perturbation, Eq. (|79p with 
5(r) — cos[ujt]. The parameters of the waiting time distri- 
bution, Eq. are a = 1/2, and = 1/2. The strength 
of the perturbation is A = 0.1. The circles correspond to an 
average over 10"^ realizations (see text). In (a) the parame- 
ters are Q = 3, uj — 1. The dotted lines are proportional to 
±Po{r). In (b), n = Lj = 1. The time r and the parameters 
are measured in arbitrary units (a.u.). 



denotes the average over the single realizations. Be- 
tween two consecutive events [action of £{t)], occurring 
at times Ti_i and t,;, the stochastic evolution is given by 
Sstir) = cos[Q.{t - Ti_i)]5'st(ri_i), where r G (n, r^.i), 
while aX T = Ti, we apply the disruptive transforma- 
tion Sst{T) A^(ri). The statistics of the time intervals 
(r^ — Ti_i) is given by the distribution w{t). In the result 
illustrated by the following figures, we use Eq. (|40p. 

In Fig. 3 we show both the solution of Eq. ([50]) and the 
average on the realizations of the stochastic simulation. 
The function /(r,0) is given by Eq. (|4ip . The external 
perturbation is ^(r) = cos[ti;T]. We note that, in general, 
in the long-time regime the system becomes insensitive to 
the external perturbation [Fig. 3a]. In fact, in the Marko- 
vian case, or when lim,-^ /(r, 0) > 0, the asymptotic 
behavior of the mean value lim7._»oo Sz(t) is given by an 
oscillatory function. In contrast, here lim,-^oo Sz{t) — 0. 
The decay to this asymptotic value is given by a power 
law function. The oscillation amplitude is proportional 
to the survival probability, Eq. which here can be 
written as Po{t) = e-y.^[A\,^t]evic[Aii2t^^'^]- In the time 
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asymptotic regime it behaves as Po{t) « 1/(^1/2%/*) 
[4J]. Only when uj = Q, the asymptotic behavior is 
given by an undamped oscillatory function [Fig. 3b]. 
The presence of an undamped asymptotic contribution 
follows straightforwardly from Eq. (|80p after expand- 
ing the involved trigonometric functions and using that 
Po(u)/(u,0) = 1/u — Po{u). The non vanishing contri- 
bution is [1 — Po{t)] cos[ljt]/2. Thus, the convergence to 
the maximal amplitude oscillation (one half) also follows 
a power law behavior. This effect is seen in Fig. 3b. 



2. Perturbing the times of events occurrence 

In Refs. [H, [3^, the response of a classical two- 
level system was analyzed by assuming that the exter- 
nal perturbation affects the times of event occurrence. 
The prediction generated by this assumption has been 
recently confirmed by experimental results on liquid crys- 
tals [4l| . Here, that assumption corresponds to assuming 
that the times of the superoperator £'s action are slightly 
changed by the external perturbation. Below, we discuss 
a system-bath modeling where this condition applies. 

We consider a system, which may also have its own 
(Markovian) dissipative dynamics, and whose interaction 
with a complex bath only occurs when the environment 
undergoes a structural change, implying the application 
of £ over the system density matrix. The changes be- 
tween the different structures of the bath are described 
by a complex landscape. The escape over a single well is 
described by the standard Kramers theory. By turning 
on an external perturbation, the height of a given well 
is written as V{t) = Vq -I- Xo^{t). By assuming an adia- 
batic regime, the survival probability associated to each 
well evolves as 



dT 



Po{r\t)~-j,{T + t)Po{T\t). 



(81) 



Here, Po{T\t) defines the conditional probability that 
no event (structural change) occurs in the time in- 
terval {t, t + t) given that the last event occurred at 
time t. The time dependent rate is written as 7o(i) = 
jQexp[—V{t)/D]. The coefficient D involves the temper- 
ature and extra parameters describing the well shape. 
The solution of Eq. ^ reads Po{T\t) ~ exp[- /J" 7o(r'-H 
t)dT']. To first order in Aq, it follows 



Po{T\t)~exp{ 



■ 7 



T + A 



'^\{t')dt']}. (82) 



where 7 = 70 exp[— Vq/I?] is the Kramers rate and the 
dimensionless strength parameter reads A = —Xo/D. If 
one assumes a statistical distribution of rate 7 (due to a 
random 7q or Vq), the survival probability Po{T\t) must 
be written as a statistical superposition of exponential 
functions, from which arbitrary decay behaviors can be 
recovered. Therefore, Eq. ([5^ can be extended to non- 
exponential survival probabilities. 



In conclusion, the external perturbation shifts the time 
of event occurrence, and this property, on the basis of the 
earlier arguments corresponds to the assumption 



Po{T\t)^Po 



AO 



t+T 



dt'({t') 



(83) 



The function ^{t') defines the time dependence of the 
perturbation. The superoperator O takes into account a 
dependence of the time shift on the system state. 
To first order in A, the survival probability reads 



PoiT\t) ^ Po{t) ~ XOw{t) 



t + T 



dt'£,{t'). 



(84) 



Although, in principle, the external perturbation breaks 
the renewal character of the process, under the assump- 
tion of weak perturbation (small A), it is legitimate to 
define the conditional waiting-time distribution density 
w{T\t) — ~{d/dT)Po{T\t), thereby getting 



w{T\t)Kw{T) + XO-^ 

OT 



w{r) i'^\t'm 



(85) 



Notice that under this approximation, independently 
of the O structure, the normalization condition 
dTw{T\t) = 1 is preserved. 
The earlier assumptions do not affect Eq. (|67|). which 
is made to remain valid by defining 0(t2, t^) as 

0(r2, ra) = O' exp[(r2 - T3)Cs]5w{t2\t3), (86) 
where O' = SO and 
d 

9t2 



5w{t2\t:>,) 



w(t2 - T3) 



dt'^{t') 



(87) 



Nevertheless, in this case, the operator Vo{t—ti) appear- 
ing in Eq. (j67p also gives a contribution to first order in 
A, as clearly shown by Eq. 

By assuming the initial condition Eq. (|75p . the condi- 
tion of Eq. ([7S)l , doing the same calculus for the deriva- 
tion of Eq. (j7ip . after rearranging the time integrals and 
some algebra, the operator expectation values read 



A(r) = Aoo + A / dT'x^o(^,^')e(r'), 
Jo 

with Aoo = Trs[Ap|^] and the response function being 
Xloir^) = w[t~t\t')Tts{A[£,0]p^} (89) 

dT"w{T" ^t',t') 

' t' 

where ■w{t — t',t') is defined by Eq. ((24l) . Thus, in this 
case the response function is also proportional to the cor- 
relation between the operator A and the external pertur- 
bation O. 
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When O'pf = or Cps" = 0, the integral contri- 
bution to Eq. (ISni) vanishes. Then Eq. ((521) yields 



A(t) = Aoo + AA'^ / dr'wir - r', t')C(t'), (90) 
Jo 

where A'^,= TTs{A[£,0]p'g'}. With the earlier argu- 
ments in mind, we state that this kind of response func- 
tion is generated whenever the underlying dynamics can 
be modeled as an escape process from a well through 
a time dependent barrier, or in general when the pre- 
scription Eq. ([55]) applies. This makes it possible to 
use for u;(t) any form, and not necessarily, the inverse 
power law form of Refs. [13, [13] ■ The example dis- 
cussed in these papers is a symmetrical two-level sys- 
tem, with £[p] = T.a,b=±^ahpalj2 = 1/2, pf = 1/2, 
C[/o] - Y.a=± ao-aapcrt^, £5 = 0, and A ^cr^. In this 
case, Xao('''' ''"') follows straightforwardly from the first 
order contribution associated to Eq. ([M]) . On the other 
hand, while Eqs. ([75)1 and ([55]) define the system response 
in terms of operator correlations, they do not involve in 
general the derivative of an operator correlation [s^l ■ 

B. Perturbing the unitary dynamics 

Now we consider the case where the external perturba- 
tion affects the unitary dynamics acting in the time in- 
terval between the occurrence of two consecutive events. 
Then, we write 

Cs{T)^Cs + \C,Ar)- (91) 

In contrast to the previous case, here we show that when 
the dynamics strongly depart from the Markovian case, 
it is not possible to generate a linear response theory. 

The perturbed dynamics, Eq. (|65p . are derived from 
the master equation defining the density matrix evolu- 
tion, Eq. with t = 0. To first order in A, we get 

'^pM^Cs{t)ps{t) + £dr'K{r - r')CU{r, r')ps{r'). 

Here, U{t' ,t") is propagator associated to the time 
dependent Liouvillian superoperator Eq. (j9ip . i.e., 
U{t', t") = exp[£s(r' - r")] + \U^^\t' ,t") + ■■■ , where 
U^^\t' ,t") is the first order contribution. 

The first-order contribution in Eq. (|65p . by assuming 
the stationary initial condition Eq. (|75p. reads 

p«(r) = J\r'g{r-T'){c^^,{T')pf (92) 

+ £ dT"K{T' - T")aA^^^ (r', t")p^ }. 

When calculating the operator expectation values, the 
first line recovers the standard Kubo response theory [sij . 
On the other hand, the second line shows that it is im- 
possible to generate a first-order perturbation. In fact. 



the validity of this contribution relies on approximating 
the difference between the perturbed and unperturbed 
propagator, U{t' , t") — exp[£s('''' — t")], by the first or- 
der contribution XU'^^^ {t' ,t"). Nevertheless, if the kernel 
K{t' — t") correlates distant times (t' and r"), evidently 
U{t' , t") cannot be approximated to first order in A. Only 
when the non-Markovian dynamics slightly depart from 
the Markovian condition [K{t' — r") = 75(r' — r")], the 
perturbed dynamics can be approximated to first order 
in the perturbation. The same conclusion follows by an- 
alyzing the perturbed stochastic trajectories. 



V. SUMMARY AND CONCLUSIONS 

In this paper we have shown that non-standard non- 
stationary statistical effects can arise in the context of 
CP open quantum system dynamics. The results rely on 
modelling the system dynamics through a renewal ap- 
proach, where the density matrix follows after averaging 
a set of realizations which mimic the interaction with a 
non-Markovian environment. The realizations are char- 
acterized by disruptive abrupt events, producing changes 
described by the application of a CP superoperator. The 
time distance between the occurrence of two consecutive 
collisional events is drawn from a non-Poisson waiting- 
time distribution density w(r). In the time intervals be- 
tween two consecutive collisional events, the system's 
time evolution is described by a unitary prescription. 
Both the CP superoperator and the waiting-time distri- 
bution density take into account the interaction of the 
system with the environment. 

As a significant advance compared to the earlier work, 
here we analyzed the non-Markovian system dynamics by 
introducing a system preparation at an arbitrary time 
and studied the ensuing evolution. The preparation 
erases the dependence of the evolution on the previous 
history of the system. Nevertheless, it does not erase the 
memory of the universe, i.e., the system-environment ar- 
rangement. In fact, the master evolution after prepara- 
tion depends explicitly on the time preparation [Eqs. (|18p 
and (|30p or Eq. pTjl ] . When the preparation time is done 
at arbitrary long times, the ensuing density matrix evo- 
lution may or not converge to an asymptotic structure, 
the last situation defining the non-stationary case. It 
arises when the average time between events is diver- 
gent. When there exists an asymptotic stationary evolu- 
tion, we showed that it may significantly depart from the 
evolution ensuing preparation at the initial time, i.e., in 
general the stationary evolution may develop stronger or 
weaker non-Markovian effects than the evolution ensuing 
the preparation at the initial time. 

The possibility of extending the regression hypothe- 
sis to the evolutions arising from the renewal approach 
was also explored. We showed that when the unitary 
dynamics commutes with the event superoperator, the 
non-Markovian evolution of expectation and correlation 
operators are exactly the same [Eqs. (|55l) and ([TO]) ]. This 
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result is valid even in the presence of non-stationary ef- 
fects. When the commutation condition is not satisfied, 
the regression hypothesis remains vahd [Eqs. ([63| and 
([M)) ] up to first order in the perturbation [Eq. (I62p ]. 

The non-stationary character of the evolution was also 
analyzed through the response of the system to an ex- 
ternal weak perturbation. When the external field mod- 
ifies the dissipative dynamics, the response function as- 
sociated to the mean value of a given operator can be 
written as a function of the correlation between the op- 
erator and the external perturbation. Different response 
functions [Eqs. ([75]) and ([55]) ] are generated depending 
on whether we make the perturbation modify the super- 
operator structure [Eq. (f68|) ] without affecting the oc- 
currence time of the collisional events or we make per- 
turbation affect the time occurrence of the disruptive 
events [Eq. ((83| ]. As in the classical counterparts, we 
have shown that in the presence of non-stationary dy- 
namics, the response of the system may die out in the 
time asymptotic regime. We also concluded that when 
the external perturbation modifies the unitary dynamics 
between events, the linear response theory is incompati- 
ble with the presence of strong memory effects. 

The equations that express the previous results are also 
valid for classical systems. In fact, all quantum proper- 
ties disappear if one disregard the unitary contributions, 
consider diagonal density matrixes, and take superoper- 
ators that do not break that condition. With respect to 



previous analysis [H, [2^, [H, [H, [H, [H, [1^ , the present 
results do not rely on a specific form of the waiting-time 
distribution (like the inverse power law forms) neither 
rely on a classical two-level system modeling. Thus, the 
formalism applies even when the events are defined by 
differential (Fokker-Planck) operators. 

The results found in this paper may have direct ex- 
perimental implications. In fact, the renewal dynamics 
arise trivially in the context of (non-Markovian) quan- 
tum kicked systems [soj. The main conclusions arrived 
at with our analysis may also apply to quantum systems 
coupled to complex reservoirs generating decay behav- 
iors without a characteristic time scale [151, [2l|. While 
a full quantum microscopic derivation of the present re- 
sults is an open problem, our analysis demonstrates that 
a reach kind of behaviors may arise when dealing with 
non-Markovian non-stationary quantum evolutions. Our 
contribution is a consistent attempt to model those issues 
in the context of CP open quantum system dynamics. 
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